function b_AFNA=Clogit_AFNA(Y_AFNA,X_AFNA,Z_AFNA)

choice=1:13;
ns=numel(choice);

Y=sortrows(Y_AFNA,[1,2]);
Y=Y(:,3:end);

X=sortrows(X_AFNA,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_AFNA=sortrows(Z_AFNA,[1,2,3]);
Z=zeros(size(Y,1),size(Z_AFNA,2),ns);
for i=4:size(Z_AFNA,2)
    Z(:,i,:)=reshape(Z_AFNA(:,i),ns,[])';
end
Z(:,1:3,:)=[];

%{
% X coefficients: constant bdtot nprofit perc_mcr hhi totpop65  
b(:,2) = [-5.866; 0.000754; -1.206; 1.748; -0.516; 0.0144];
b(:,3) = [-8.014; 0.00104; 0.493; -0.189; 0.708; 0.225];
b(:,4) = [-6.345; -0.0118; -0.727; 1.976; -0.323; 0.0959 ];
b(:,5) = [-3.321; -0.0124; -0.394; 0.540; 0.515; -0.181 ];
b(:,6) = [-13.54; 0.000318; 1.350; -0.419; 1.209; 0.551 ];
b(:,7) = [-8.657; -0.0000928; 3.153; -0.442; -0.297; 0.0860 ];
b(:,8) = [-12.86; 0.00113; 2.111; -0.556; 2.096; 0.415 ];
b(:,9) = [-4.742; -0.00298; -1.673; -1.117; 1.376; 0.0365 ];
b(:,10) = [-7.240; 0.000858; 0.131; -0.761; 0.421; 0.201 ];
b(:,11) = [-7.661; -0.000438; 0.331; -0.908; 0.654; 0.225 ];
b(:,12) = [-3.975; 0.000417; -0.320; -0.480; -0.291; -0.00772 ];
b(:,13) = [-7.754; 0.000561; 0.560; -2.204; 0.881; 0.230];
%}


% ### UPDATED on 06/08/2020: only keep bdtot, HHI, and totpop65
% X coefficients: constant bdtot hhi totpop65  
b(:,2) = [-5.866; 0.000754;  -0.516; 0.0144];
b(:,3) = [-8.014; 0.00104; 0.708; 0.225];
b(:,4) = [-6.345; -0.0118; -0.323; 0.0959 ];
b(:,5) = [-3.321; -0.0124;  0.515; -0.181 ];
b(:,6) = [-13.54; 0.000318;  1.209; 0.551 ];
b(:,7) = [-8.657; -0.0000928;  -0.297; 0.0860 ];
b(:,8) = [-12.86; 0.00113;  2.096; 0.415 ];
b(:,9) = [-4.742; -0.00298; 1.376; 0.0365 ];
b(:,10) = [-7.240; 0.000858;  0.421; 0.201 ];
b(:,11) = [-7.661; -0.000438;  0.654; 0.225 ];
b(:,12) = [-3.975; 0.000417;  -0.291; -0.00772 ];
b(:,13) = [-7.754; 0.000561; 0.881; 0.230];


% Z coefficients: vsysh, MKTabove25%, MKTabove50%, MKTabove75%
bz=[2.788; 0.284; 0.535; 1.011];

bOri=[reshape(b(:,2:end),[],1);bz];

%bAns=zeros(size(bOri));

options=optimset('MaxFunEvals',200000000000,'MaxIter',1500000000000,'TolX',1e-16,'Tolfun',1e-16,'GradObj','on','DerivativeCheck','off','FinDiffType','central');
startval = bOri;
baseAlt=1;
b_AFNA = fminunc('clogit',startval,options,[],Y,X,Z,baseAlt);













